Revision: T-cell activation
3/13/2020
1 Introduction
Reviewer #2 requested the following analysis for the T-cell activation and culturing experiment:
- Characterization of the transcriptomic changes associated with exposure to anti-CD3 antibody.
- Compositional analysis of cell types with and without anti-CD3 antibody treatment.
2 Pre-processing
2.1 Package loading
library(scater)
library(Seurat)
library(ggpubr)
library(purrr)
library(kBET)
library(viridis)
library(ggmap)
library(cowplot)
library(tidyverse)2.2 Source script with functions
source("bin/utils.R")2.3 Load data
t_act_male1 <- readRDS("../3-T_cell_activation/results/R_objects/t_act_Seurat_male1.rds")
t_act_rep2_l <- readRDS("../3-T_cell_activation/results/R_objects/t_act_Seurat_list_reannotated_rep2.rds")2.4 Merge seurat objects and homogenize annotations
# Donor 1
t_act_male1$`0M`$cell_type <- factor(
t_act_male1$`0M`$cell_type,
levels = c("CD4 T", "Cytotoxic", "Monocyte", "B", "FCGR3A Monocyte", "Dendritic Cell")
)
levels(t_act_male1$`0M`$cell_type) <- c("CD4 T-cell", "Cytotoxic", "CD14 Monocyte", "B-cell",
"FCGR3A Monocyte", "Dendritic Cell")
t_act_male1$`2M`$cell_type <- factor(
t_act_male1$`2M`$cell_type,
levels = c("CD4 T", "Cycling", "Cytotoxic", "B")
)
levels(t_act_male1$`2M`$cell_type) <- c("Activated CD4 T-cell", "Cycling CD4 T-cell",
"Cytotoxic", "B-cell")
t_act_male1$`0M`$is_cultured <- rep("uncultured", ncol(t_act_male1$`0M`))
t_act_male1$`2M`$is_cultured <- rep("cultured", ncol(t_act_male1$`2M`))
t_act_male1 <- merge(
x = t_act_male1$`0M`,
y = t_act_male1$`2M`,
add.cell.ids = c("uncultured", "cultured")
)
t_act_male1$donor <- rep("donor1", ncol(t_act_male1))
# Donor 2
t_act_rep2_l$`0_rep2_F2`$cell_type <- factor(
t_act_rep2_l$`0_rep2_F2`$cell_type,
levels = c("CD4 T-cell", "CD8 T-cell", "NK", "B-cell", "Monocyte", "FCGR3A Monocyte", "Unknown")
)
levels(t_act_rep2_l$`0_rep2_F2`$cell_type) <- c("CD4 T-cell", "CD8 T-cell", "NK",
"B-cell", "CD14 Monocyte", "FCGR3A Monocyte",
"Unknown")
t_act_rep2_l$`1_rep2_F2`$cell_type <- factor(
t_act_rep2_l$`1_rep2_F2`$cell_type,
levels = c("Activated CD4 T-cell", "Cycling CD4 T-cell", "CD8 T-cell", "NK", "B-cell", "Unknown")
)
t_act_rep2_l$`0_rep2_F2`$is_cultured <- rep("uncultured", ncol(t_act_rep2_l$`0_rep2_F2`))
t_act_rep2_l$`1_rep2_F2`$is_cultured <- rep("cultured", ncol(t_act_rep2_l$`1_rep2_F2`))
t_act_female2 <- merge(
x = t_act_rep2_l$`0_rep2_F2`,
y = t_act_rep2_l$`1_rep2_F2`,
add.cell.ids = c("uncultured", "cultured")
)
t_act_female2$donor <- rep("donor2", ncol(t_act_female2))
# Donor 3
t_act_rep2_l$`0_rep2_F3`$cell_type <- factor(
t_act_rep2_l$`0_rep2_F3`$cell_type,
levels = c("CD4 T-cell", "CD8 T-cell", "NK", "B-cell", "Monocyte", "FCGR3A Monocyte")
)
levels(t_act_rep2_l$`0_rep2_F3`$cell_type) <- c("CD4 T-cell", "CD8 T-cell", "NK",
"B-cell", "CD14 Monocyte", "FCGR3A Monocyte")
t_act_rep2_l$`1_rep2_F3`$cell_type <- factor(
t_act_rep2_l$`1_rep2_F3`$cell_type,
levels = c("Activated CD4 T-cell", "Cycling CD4 T-cell", "CD8 T-cell", "NK", "B-cell", "Unknown")
)
t_act_rep2_l$`0_rep2_F3`$is_cultured <- rep("uncultured", ncol(t_act_rep2_l$`0_rep2_F3`))
t_act_rep2_l$`1_rep2_F3`$is_cultured <- rep("cultured", ncol(t_act_rep2_l$`1_rep2_F3`))
t_act_female3 <- merge(
x = t_act_rep2_l$`0_rep2_F3`,
y = t_act_rep2_l$`1_rep2_F3`,
add.cell.ids = c("uncultured", "cultured")
)
t_act_female3$donor <- rep("donor3", ncol(t_act_female3))
# Create list
t_act_l <- list(t_act_male1, t_act_female2, t_act_female3)
names(t_act_l) <- c("donor1", "donor2", "donor3")
t_act_l$donor1$cell_type <- factor(
t_act_l$donor1$cell_type,
levels = c("CD4 T-cell", "Activated CD4 T-cell", "Cycling CD4 T-cell", "Cytotoxic",
"B-cell", "CD14 Monocyte", "FCGR3A Monocyte", "Dendritic Cell")
)
levels_type_female <- c("CD4 T-cell", "Activated CD4 T-cell", "Cycling CD4 T-cell",
"CD8 T-cell", "NK", "B-cell", "CD14 Monocyte", "FCGR3A Monocyte",
"Unknown")
t_act_l$donor2$cell_type <- factor(
t_act_l$donor2$cell_type,
levels = levels_type_female
)
t_act_l$donor3$cell_type <- factor(
t_act_l$donor3$cell_type,
levels = levels_type_female
)
t_act_sub_l <- purrr::map(t_act_l, function(seurat) {
Idents(seurat) <- "time"
seurat_sub <- subset(seurat, idents = "0h")
seurat_sub <- seurat_sub %>%
FindVariableFeatures() %>%
ScaleData() %>%
RunPCA() %>%
RunTSNE(reduction = "pca", dims = 1:30) %>%
RunUMAP(reduction = "pca", dims = 1:30)
Idents(seurat_sub) <- "cell_type"
seurat_sub
})3 tSNE uncultured/cultured
donors <- c("Donor 1", "Donor 2", "Donor 3")
tsnes_is_cultured <- purrr::map2(t_act_sub_l, donors, function(seurat, titl) {
Idents(seurat) <- "is_cultured"
tsne <- DimPlot(
seurat,
reduction = "tsne",
pt.size = 0.5,
cols = c("#3374A1", "#E1812C")
)
tsne +
labs(title = titl, x = "tSNE1", y = "tSNE2") +
scale_color_manual("", values = c("#3374A1", "#E1812C"),
labels = c("Original", "Cultured")) +
theme(plot.title = element_text(size = 13, hjust = 0.5),
axis.ticks = element_blank(),
axis.text = element_blank(),
axis.title = element_blank(),
axis.line = element_blank())
})
tsnes_is_cultured_arr <- ggarrange(
plotlist = tsnes_is_cultured,
ncol = 3,
common.legend = TRUE
)
tsnes_is_cultured_arr
# Get legend
p <- tsnes_is_cultured$donor1 + theme(legend.position = "bottom")
leg <- as_ggplot(get_legend(p))
date <- Sys.Date()
# ggsave(
# filename = str_c("../doc/figures/legends/", date, "_", "is_cultured", ".pdf"),
# plot = leg,
# width = 16,
# height = 5,
# units = "cm"
# )4 tSNE cell types
all_cell_types <- c("CD4 T-cell", "Activated CD4 T-cell", "Cycling CD4 T-cell", "Cytotoxic",
"CD8 T-cell", "NK", "B-cell", "CD14 Monocyte", "FCGR3A Monocyte",
"Dendritic Cell", "Unknown")
all_cell_types <- factor(all_cell_types, levels = all_cell_types)
palette <- c("#81324c", "#b82e57", "#e04d74", "#3a2edb", "#752bbf",
"#c03fe4", "#bbaa2a", "#71bdd0", "green4", "hotpink2",
"gray50")
names(palette) <- all_cell_types
tsnes_cell_type <- purrr::map(t_act_sub_l, function(seurat) {
tsne <- DimPlot(
seurat,
reduction = "tsne",
pt.size = 0.5,
cols = palette[levels(seurat$cell_type)]
)
tsne +
labs(x = "tSNE1", y = "tSNE2") +
theme(axis.ticks = element_blank(),
axis.text = element_blank(),
axis.title = element_blank(),
axis.line = element_blank())
})
tsnes_cell_type_arr <- ggarrange(
plotlist = tsnes_cell_type,
ncol = 3,
common.legend = TRUE
)
tsnes_cell_type_arr
5 Cell cycle score
Garcia-Sousa I, et al. showed that the major biological process activated upon treatment with an anti-CD3 antibody is cell cycle. Let us score each cell with a cell cycle transcriptomic signature:
tsnes_s_score <- purrr::map(t_act_sub_l, function(seurat) {
tsne <- FeaturePlot(
seurat,
reduction = "tsne",
feature = "S.Score",
pt.size = 0.5,
cols = viridis(5)
)
tsne +
labs(x = "tSNE1", y = "tSNE2") +
theme(axis.ticks = element_blank(),
axis.text = element_blank(),
axis.title = element_blank(),
axis.line = element_blank(),
plot.title = element_blank())
})
tsnes_s_score_arr <- ggarrange(plotlist = tsnes_s_score, ncol = 3)
tsnes_s_score_arr
# Get legend
p <- tsnes_s_score$donor1
leg <- as_ggplot(get_legend(p))
# ggsave(
# filename = str_c("../doc/figures/legends/", date, "_", "s_phase_score", ".pdf"),
# plot = leg,
# width = 16,
# height = 5,
# units = "cm"
# )6 Compositional analysis
donors <- c("Donor 1", "Donor 2", "Donor 3")
compositional_analysis_gg <- purrr::map(t_act_sub_l, function(seurat) {
df <- seurat@meta.data %>%
dplyr::select("is_cultured", "cell_type") %>%
group_by(is_cultured, cell_type) %>%
summarise(n_cells = n()) %>%
ungroup() %>%
group_by(is_cultured) %>%
mutate(pct_cells = n_cells / sum(n_cells) * 100)
df$is_cultured <- factor(df$is_cultured, levels = c("uncultured", "cultured"))
p <- df %>%
ggplot(aes(is_cultured, pct_cells, fill = cell_type)) +
geom_col() +
labs(x = "", y = "Percentage of cells (%)", fill = "") +
scale_x_discrete(labels = c("Original", "Cultured")) +
scale_fill_manual(values = palette[levels(seurat$cell_type)]) +
theme_classic() +
theme(plot.title = element_text(size = 13, hjust = 0.5),
axis.title.y = element_text(size = 12),
axis.text.x = element_text(size = 11, color = "black"))
p
})
compositional_analysis_arr <- ggarrange(
plotlist = compositional_analysis_gg,
ncol = 3,
common.legend = TRUE
)
compositional_analysis_arr
# Get legend
p1 <- compositional_analysis_gg$donor1 +
theme(legend.position = "bottom")
p2 <- compositional_analysis_gg$donor2 +
theme(legend.position = "bottom")
leg1 <- as_ggplot(get_legend(p1))
leg2 <- as_ggplot(get_legend(p2))
leg_list <- list(cell_types1 = leg1, cell_types2 = leg2)
date <- Sys.Date()
# walk(names(leg_list), function(leg) {
# ggsave(
# filename = str_c("../doc/figures/legends/", date, "_", leg, ".pdf"),
# plot = leg_list[[leg]],
# width = 16,
# height = 5,
# units = "cm"
# )
# })7 Arrange figure
fig <- plot_grid(tsnes_is_cultured_arr, NULL, tsnes_cell_type_arr, NULL, tsnes_s_score_arr,
NULL, compositional_analysis_arr, nrow = 7, ncol = 1,
rel_heights = c(1, 0.05, 1, 0.1, 1, 0.05, 1))
# fig
# ggsave(filename = "../doc/figures/R/suppZZ.pdf", plot = fig, width = 18.5, height = 27.5, units = "cm")8 Differential expression analysis (cultured VS uncultured)
# Merge
t_act <- merge(
x = t_act_l$donor1,
y = c(t_act_l$donor2, t_act_l$donor3)
)
# Process
t_act <- t_act %>%
FindVariableFeatures() %>%
ScaleData() %>%
RunPCA() %>%
RunTSNE(reduction = "pca", dims = 1:30) %>%
RunUMAP(reduction = "pca", dims = 1:30)
Idents(t_act) <- "donor"
DimPlot(t_act)
t_act$cell_type[t_act$cell_type %in% c("NK", "CD8 T-cell")] <- "Cytotoxic"
Idents(t_act) <- "cell_type"
t_act <- subset(
t_act,
idents = c("CD4 T-cell", "Activated CD4 T-cell", "Cytotoxic", "B-cell")
)
# DEA
t_act$cell_type[t_act$cell_type == "Activated CD4 T-cell"] <- "CD4 T-cell"
Idents(t_act) <- "cell_type"
types <- c("CD4 T-cell", "Cytotoxic", "B-cell")
t_act_types <- SplitObject(t_act, split.by = "cell_type")
dea_l <- purrr::map2(t_act_types, types, function(seurat, type) {
Idents(seurat) <- "is_cultured"
df <- FindMarkers(
seurat,
ident.1 = "cultured",
ident.2 = "uncultured",
test.use = "wilcox",
logfc.threshold = 0
)
df
})
dea_df_l <- purrr::map2(dea_l, t_act_types, function(df, seurat) {
mat <- as.matrix(seurat[["RNA"]]@data[rownames(df), ])
average_expression <- rowMeans(mat)
log2_fc <- apply(mat, 1, function(x) {
mean_uncultured <- mean(x[seurat$is_cultured == "uncultured"]) + 0.05
mean_cultured <- mean(x[seurat$is_cultured == "cultured"]) + 0.05
log2(mean_cultured / mean_uncultured)
})
df <- df %>%
rownames_to_column(var = "gene") %>%
dplyr::select("gene", "p_val", "p_val_adj") %>%
dplyr::mutate(average_expression = average_expression, log2_fc = log2_fc) %>%
dplyr::arrange(p_val_adj) %>%
dplyr::filter(p_val_adj < 0.001)
df <- df[, c("gene", "average_expression", "log2_fc", "p_val", "p_val_adj")]
df
})
DT::datatable(dea_df_l$`CD4 T-cell`)DT::datatable(dea_df_l$Cytotoxic)DT::datatable(dea_df_l$`B-cell`)# openxlsx::write.xlsx(dea_df_l, file = "results/tables/dea_culturedVSuncultured.xlsx")9 Session Information
sessionInfo()## R version 3.6.1 (2019-07-05)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] forcats_0.5.0 stringr_1.4.0 dplyr_0.8.4 readr_1.3.1 tidyr_1.0.2 tibble_2.1.3 tidyverse_1.3.0 cowplot_1.0.0 ggmap_3.0.0 viridis_0.5.1 viridisLite_0.3.0 kBET_0.99.6 purrr_0.3.3 ggpubr_0.2.5 magrittr_1.5 Seurat_3.1.4 scater_1.14.6 ggplot2_3.3.0 SingleCellExperiment_1.8.0 SummarizedExperiment_1.16.1 DelayedArray_0.12.2 BiocParallel_1.20.1 matrixStats_0.55.0 Biobase_2.46.0 GenomicRanges_1.38.0 GenomeInfoDb_1.22.0 IRanges_2.20.2 S4Vectors_0.24.3 BiocGenerics_0.32.0 BiocStyle_2.14.4
##
## loaded via a namespace (and not attached):
## [1] reticulate_1.14 tidyselect_1.0.0 htmlwidgets_1.5.1 grid_3.6.1 Rtsne_0.15 munsell_0.5.0 codetools_0.2-16 mutoss_0.1-12 ica_1.0-2 DT_0.12 future_1.16.0 withr_2.1.2 colorspace_1.4-1 knitr_1.28 rstudioapi_0.11 ROCR_1.0-7 ggsignif_0.6.0 gbRd_0.4-11 listenv_0.8.0 Rdpack_0.11-1 labeling_0.3 RgoogleMaps_1.4.5.3 GenomeInfoDbData_1.2.2 mnormt_1.5-6 farver_2.0.3 vctrs_0.2.3 generics_0.0.2 TH.data_1.0-10 xfun_0.12 R6_2.4.1 ggbeeswarm_0.6.0 rsvd_1.0.3 bitops_1.0-6 assertthat_0.2.1 promises_1.1.0 scales_1.1.0 multcomp_1.4-12 beeswarm_0.2.3 gtable_0.3.0 npsurv_0.4-0 globals_0.12.5 sandwich_2.5-1 rlang_0.4.5 splines_3.6.1 lazyeval_0.2.2 broom_0.5.5 BiocManager_1.30.10
## [48] yaml_2.2.1 reshape2_1.4.3 modelr_0.1.6 crosstalk_1.0.0 backports_1.1.5 httpuv_1.5.2 tools_3.6.1 bookdown_0.18 gplots_3.0.3 RColorBrewer_1.1-2 ggridges_0.5.2 TFisher_0.2.0 Rcpp_1.0.3 plyr_1.8.6 zlibbioc_1.32.0 RCurl_1.98-1.1 pbapply_1.4-2 zoo_1.8-7 haven_2.2.0 ggrepel_0.8.1 cluster_2.1.0 fs_1.3.2 data.table_1.12.8 RSpectra_0.16-0 magick_2.3 lmtest_0.9-37 reprex_0.3.0 RANN_2.6.1 mvtnorm_1.1-0 fitdistrplus_1.0-14 xtable_1.8-4 mime_0.9 hms_0.5.3 patchwork_1.0.0 lsei_1.2-0 evaluate_0.14 jpeg_0.1-8.1 readxl_1.3.1 gridExtra_2.3 compiler_3.6.1 KernSmooth_2.23-16 crayon_1.3.4 htmltools_0.4.0 later_1.0.0 RcppParallel_4.4.4 lubridate_1.7.4 DBI_1.1.0
## [95] dbplyr_1.4.2 MASS_7.3-51.5 Matrix_1.2-18 cli_2.0.2 gdata_2.18.0 metap_1.3 igraph_1.2.4.2 pkgconfig_2.0.3 sn_1.5-5 numDeriv_2016.8-1.1 sp_1.4-1 plotly_4.9.2 xml2_1.2.2 vipor_0.4.5 multtest_2.42.0 XVector_0.26.0 bibtex_0.4.2.2 rvest_0.3.5 digest_0.6.25 sctransform_0.2.1 RcppAnnoy_0.0.15 tsne_0.1-3 rmarkdown_2.1 cellranger_1.1.0 leiden_0.3.3 uwot_0.1.5 DelayedMatrixStats_1.8.0 shiny_1.4.0 gtools_3.8.1 rjson_0.2.20 lifecycle_0.1.0 nlme_3.1-145 jsonlite_1.6.1 BiocNeighbors_1.4.2 fansi_0.4.1 pillar_1.4.3 lattice_0.20-40 fastmap_1.0.1 httr_1.4.1 plotrix_3.7-7 survival_3.1-8 glue_1.3.1 FNN_1.1.3 png_0.1-7 stringi_1.4.6 BiocSingular_1.2.2 caTools_1.18.0
## [142] irlba_2.3.3 future.apply_1.4.0 ape_5.3